// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Copyright (C) December 9, 2009 MARIN - All rights reserved - http://www.marin.nl
// -----------------------------------------------------------------------------
// Program  :
// Module   :
// File     : ShipPulse.cpp
// -----------------------------------------------------------------------------

#include <new>

#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "defines.h"
#include <mml/Vector.h>
#include "ShipPulse.h"
#include "Timer.h"
#include "assertplus.h"

// for testing purposes, these are global variables
int safe = 0;
int jj = 0;
int shipPulseCount = 0;
double shipPulseTime = 0;

// -----------------------------------------------------------------------------
// function determineBoundingBoxCells
// -----------------------------------------------------------------------------
/// determines cells corresponding to bounding box vertices
// -----------------------------------------------------------------------------
//static
void determineBoundingBoxCells(double a_xmin, ///< minimum x coordinate of panel bounding box
                               double a_ymin, ///< minimum y coordinate of panel bounding box
                               double a_xmax, ///< maximum x coordinate of panel bounding box
                               double a_ymax, ///< maximum y coordinate of panel bounding box
                               double a_dx,   ///< ew grid spacing
                               double a_dy,   ///< ns grid spacing
                               int* a_pimin,  ///< minimum ew grid point index
                               int* a_pjmin,  ///< minimum ns grid point index
                               int* a_pimax,  ///< maximum ew grid point index
                               int* a_pjmax)  ///< maximum ns grid point index
{
    a_xmin += 0.5 * a_dx;
    a_ymin += 0.5 * a_dy;
    a_xmax += 0.5 * a_dx;
    a_ymax += 0.5 * a_dy;

    double x, y;
    modf(a_xmin / a_dx, &x);
    modf(a_ymin / a_dy, &y);
    *a_pimin = int(x);
    *a_pjmin = int(y);
    modf(a_xmax / a_dx, &x);
    modf(a_ymax / a_dy, &y);
    *a_pimax = int(x);
    *a_pjmax = int(y);
#ifdef SMOOTHING  // TODO!!!
    *a_pimin = *a_pimin - 1;
    *a_pjmin = *a_pjmin - 1;
    *a_pimax = *a_pimax + 1;
    *a_pjmax = *a_pjmax + 1;
#endif
}

bool ShipPulse::intersect(const xxship::Panel& a_panel,   ///< panel to cut up
                          int a_i,                        ///< ew index of grid point
                          int a_j,                        ///< ns index of grid point
                          xxship::Panel* a_pIntersection) ///< new panel with intersection result
{
    const mml::Planed* pPlanes[4];
    int i = a_i - m_ewOffset;
    int j = a_j - m_nsOffset;
    pPlanes[0] = &(m_wplanes[i][j]);
    pPlanes[1] = &(m_splanes[i][j]);
    pPlanes[2] = &(m_eplanes[i][j]);
    pPlanes[3] = &(m_nplanes[i][j]);
    *a_pIntersection = a_panel.createReducedPanel(pPlanes);
    return a_pIntersection->getNumberOfVertices() >= 3;
}

// -----------------------------------------------------------------------------
// function ShipPulse::ShipPulse
// -----------------------------------------------------------------------------
/// constructor
/** detailed description */
// -----------------------------------------------------------------------------
ShipPulse::ShipPulse(const ShipGeometry& a_ship, double a_dx, double a_dy)
: super(a_dx, a_dy)
, m_geometry(a_ship)
{
}


// -----------------------------------------------------------------------------
// function ShipPulse::~ShipPulse()
// -----------------------------------------------------------------------------
/// destructor
/** detailed description */
// -----------------------------------------------------------------------------
ShipPulse::~ShipPulse()
{
    // implementation
}


// -----------------------------------------------------------------------------
// function ShipPulse::compute
// -----------------------------------------------------------------------------
/// computes pressure pulse for use with Gert Klopman's wave model
// -----------------------------------------------------------------------------
void ShipPulse::compute(int *a_px,   ///< x offset of resulting pulse array
                        int *a_py,   ///< y offset of resulting pulse array
                        Array2D<REAL>* a_pPulse) ///< pulsedepths
{
    double t = seconds();

    // syntactic sugar
    Array2D<REAL>& pulse = *a_pPulse;

    // determine indices of horizontal bounding box corners in waveheight array
    const xxship::Panels& earthPanels = m_geometry.getEarthPanels();
    double xmin;
    double ymin;
    double xmax;
    double ymax;
    earthPanels.getBoundingBox(&xmin, &ymin, &xmax, &ymax);

    // make potatocutter for entire bounding box
    int imin;
    int jmin;
    int imax;
    int jmax;
    determineBoundingBoxCells(xmin, ymin, xmax, ymax,
                              m_dx, m_dy,
                              &imin, &jmin, &imax, &jmax);
    makePlanes(imin, jmin, imax, jmax);

    // set dimensions of pulse array to size of bounding box
    pulse.clear();
    int dim1 = imax - imin + 1;
    int dim2 = jmax - jmin + 1;
    pulse.resize(dim1, dim2);

    // determine offsets of pulse array within waveheight array
    *a_px = imin;
    *a_py = jmin;

    // set all depth values in pulse array to 0
    pulse = 0.0;

    // loop over all panels
    int n = m_geometry.getNumberOfPanels();
#ifdef openmp3
#pragma omp parallel for
#endif
    for (int p = 0; p < n; ++p)
    {
        const xxship::Panel& panel = earthPanels.getPanel(p);
        if (panel.getMinZ() < 0) // only look at panels below 0
        {
            // necessary for bulb ships:
            // determine normal on panel, if pointing below XY plane
            // add contribution of plane, otherwise subtract it
            // (normals point outward from panel into water)
            int sign = panel.getSurfaceVector() * mml::Vector3d(0,0,-1) < 0 ? -1 : 1;

            // determine grid indices of panel's bounding box
            int iimin, iimax, jjmin, jjmax;
            double xxmin, yymin, xxmax, yymax;
            panel.getBoundingBox(&xxmin, &yymin, &xxmax, &yymax);
            determineBoundingBoxCells(xxmin, yymin,
                                      xxmax, yymax,
                                      m_dx, m_dy,
                                      &iimin, &jjmin, &iimax, &jjmax);

            if (iimin == iimax && jjmin == jjmax) // panel inside one grid cell; no chopping necessary
            {
                double depth = panel.getCentroid()(2);
                if (depth < 0)
                {
                    int ii = iimin - imin;
                    int jj = jjmin - jmin;
#ifdef openmp3
#pragma omp critical
#endif
                    { // -----------------------------------------------------------------------
                        // add contribution of panel to pulsedepth in gridcell:
                        // XY-projected-area-of-chopped-panel * depth-at-centroid-of-chopped-panel / total-area-of-cell
                        pulse[ii][jj] += sign * panel.getProjectedArea() * depth * m_invertedCellArea;
                    } // critical --------------------------------------------------------------
                }
            }
            else
            {
                // loop over all grid cells in bounding box with 'potatoCutter'
                for (int i = iimin; i <= iimax; ++i)
                {
                    for (int j = jjmin; j <= jjmax; ++j)
                    {
                        // clip panel by one 'potatoCutter' cell
                        xxship::Panel result; // receives the clipping result
                        if (intersect(panel, i, j, &result))
                        {
                            double depth = result.getCentroid()(2);
                            if (depth < 0)
                            {
                                int ii = i - imin;
                                int jj = j - jmin;
#ifdef openmp3
#pragma omp critical
#endif
                                { // -----------------------------------------------------------------------
                                    // add contribution of panel to pulsedepth in gridcell:
                                    // XY-projected-area-of-chopped-panel * depth-at-centroid-of-chopped-panel / total-area-of-cell
                                    pulse[ii][jj] += sign * result.getProjectedArea() * depth * m_invertedCellArea;
                                } // critical --------------------------------------------------------------
                            }
                        }
                        else
                        {
//                            ASSERT(false, "");
                        }
                    }
                }
            }
        }
    }
#ifdef PRINT
    PrintArray2Dcontents("", pulse, 0, dim1 - 1, 0, dim2 - 1);
#endif
    t = seconds() - t;
    shipPulseTime += t;
    shipPulseCount++;
}

void ShipPulse::makePlanes(int a_imin, int a_jmin, int a_imax, int a_jmax)
{
    int ew = a_imax - a_imin + 1;
    int ns = a_jmax - a_jmin + 1;
    m_ewOffset = a_imin;
    m_nsOffset = a_jmin;
    m_nplanes.resize(ew, ns);
    m_splanes.resize(ew, ns);
    m_eplanes.resize(ew, ns);
    m_wplanes.resize(ew, ns);
#ifdef openmp3
#pragma omp parallel for
#endif
    for (int i = 0; i < ew; ++i)
    {
        for (int j = 0; j < ns; ++j)
        {
            double w = (i + a_imin - 0.5) * m_dx;
            double s = (j + a_jmin - 0.5) * m_dy;
            double e = (i + a_imin + 0.5) * m_dx;
            double n = (j + a_jmin + 0.5) * m_dy;
            m_wplanes[i][j] = mml::Planed(mst::Vector3d(1, 0, 0), w);
            m_splanes[i][j] = mml::Planed(mst::Vector3d(0, 1, 0), s);
            m_eplanes[i][j] = mml::Planed(mst::Vector3d(-1, 0, 0), -e);
            m_nplanes[i][j] = mml::Planed(mst::Vector3d(0, -1, 0), -n);
        }
    }
}
